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ABSTRACT 

Context. The increasing interest in understanding stellar magnetic activity cycles is a strong motivation for the development of 
parameterized starspot models which can be constrained observationally. 

Aims. In this work we develop an empirical tool for the stochastic reconstruction of sunspot cycles, using the average solar properties 
as a reference. 

Methods. The synthetic sunspot cycle is compared with the sunspot data extracted from the National Geophysical Data Center, in 
particular using the Kolmogorov-Smirnov test. This tool yields synthetic spot group records, including date, area, latitude, longitude, 
rotation rate of the solar surface at the group’s latitude, and an identification number. 

Results. Comparison of the stochastic reconstructions with the daily sunspot records confirms that our empirical model is able to 
successfully reproduce the main properties of the solar sunspot cycle. As a by-product of this work, we show that the Gnevyshev- 
Waldmeier rule, which describes the spots’ area-lifetime relation, is not adequate for small groups and we propose an effective 
correction to that relation which leads to a closer agreement between the synthetic sunspot cycle and the observations. 
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1. Introduction 


In the Sun, the level of magnetic activity varies over time show¬ 
ing a periodic behaviour known as the solar cycle. The most di¬ 
rect evidence for this is the 11-year variation in the number of 
sunspots observed at the solar photosphere. As the cycle pro¬ 
ceeds, the sunspot formation zone gradually migrates towards 
the equator until the next minimum is reached and a new cycle 
begins. At that point, the polarity of the magnetic field reverses 
and returns to the original state at the end of the second 11-year 
cycle, thus completing a 22-year Hale cycle (the sunspot polarity 
law - |Hale et al.|19l9{|Hale & Nicholson|1925) . 

Evidence for the presence of activity cycles in other solar¬ 
like stars, including exoplanet hosts, has also been accumulating 
over the past years. Starspots cannot be observed directly at the 
surface of distant stars, but can be detected indirectly through 
the effects they induce, such as the strong emission at the cen¬ 
tre of the Can H and K lines. By monitoring this emission, it is 
possible to infer the rotation period of the star and the period of 


& Baliunas 1995 

Saar & Brandenburg 

1999; Cincunegui et al. 

2007; Hall et al. 

2007 ;lMetcalfe et al. 

2010,12013). Moreover, 


since spots are darker than the average stellar surface they can 
also be detected through the inspection of the photometric light 
curves of stars observed by CoRoT ((Baglin et al.|2006| and Ke¬ 
pler (Borucki et al.|[20~i~0l Koch et al.||2010|> space missions. In 


an active star with starspots crossing the visible disk, the light 
curve shows a quasi-periodic modulation which results from the 


combination of the effects of stellar rotation and magnetic ac¬ 
tivity (e.g.|Mosser et al.|2009[ Mathur et al.|2010 Garcia et al.| 


2010b| Ballot et al.|201 1| Campante et al.|2011| >. That modula¬ 

tion is more significant in periods of maximum activity, making 
the light curves a possible starspot proxy. In addition, through 
the wavelet analysis of the light curves it might be possible to 
obtain information about the time evolution of the starspots at 
the stellar surface and, in turn, derive the period of the activ¬ 
ity cycle and the rotation period of the star (e.g. Campante 2012; 


Math ur et al.|2014| |Garclaet al.|2014[ |Bravo et al.|2014| ) .Finally, 
stellar activity cycles may also be detected through their impact 
on stellar oscillations because the magnetic activity affects the 
wave propagation, inducing changes to the oscillation frequen¬ 
cies, amplitudes, and line widths, which are thus found to vary 


in phase with other sun-/starspot proxies (e.g. Woodard & Noyes 


Libbrecht & Woodard 1990; Chaplin et al.|2004[|Metcalfe 


1985 


2007| Garcia et al.|2010a| Tripathy et al.|201 1 1 . 


et al. 


Starspots simulations can be used to study how the activity- 
induced variability in the light curve and in the oscillation prop¬ 
erties depends on the characteristics of the activity cycle. These 
simulations are important in other astrophysical contexts, such 
as in the quest for exoplanets, where they can be used to design 
new strategies to reduce the signatures induced by stellar activ- 


ity on the radial velocity and transit observations (e.g. Pont et al. 

2008 

Czesla et al. 2009; 

Figueira et al.|2010, Dumusque et al. 

2011 

Oshagh et al. 12013 

2014). With the above in mind, in this 


work we develop a parameterized model aimed at reproducing 
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the main properties of the sunspot cycle that can also be applied 
to reproducing the activity cycles of other stars. In order to build 
our empirical model, we use as inputs a number of properties of 
the observed sunspot cycle. These properties are summarized in 
section |2] and their implementation in our model is discussed in 
section R] The results obtained with this tool and their compari¬ 
son with the solar data extracted from the National Geophysical 
Data Center (NOAA/NGDC) are presented in section[4] Finally, 
in section[5]we draw our main conclusions. 


2. Properties of the solar cycle 

In order to reproduce the sunspot cycle, a number of important 
observational constraints must be considered. In what follows 
we review key observational properties of the sunspot cycle that 
will be used in our model, namely: the number of sunspot groups 
and its dependence on the phase of the cycle; the sunspots’ areas 
and their relation with the sunspots’ lifetimes; the formation lat¬ 
itude of the sunspots and the width of the formation region; the 
differential rotation of the solar surface. 

The number of sunspots and sunspot groups varies over the 
solar cycle. Its evolution is asymmetric: the rising phase of the 
cycle is faster than the declining towards the next minimum. 
Different authors have used different functions to describe the 
asymmetric shape of the sunspot cycle (e.g.|Setwart & Panofsky 


N s (t) = 


ai(t - to) 3 


exp((r - to) 2 lb\) - ci 


(1) 


where to is the starting time (about four months prior to the min¬ 
imum for an average cycle ( jHathawa~y||2010f >), a\ is the ampli¬ 
tude, b\ is related to the size of the cycle, and ci is related to the 
asymmetry of the cycle. 

The asymmetric shape of the solar cycle is also evident in 
the temporal variation of the sunspot areas. Bogdan et al. ( 1988) 
were the first to notice that the accumulated umbral areas dis¬ 
tribution (derived from daily records by counting each spot as 
many times as the number of days it remains visible) can be 
described by a log-normal distribution. Since the ratio between 
the umbral area and the total area of the spot does not depend 
significantly on the spot size (e.g. Brandt et al. 1990| [Solanki 
|2003| |Vaquero et al.||2005| |Kiess et al.||2014| >7~the accumulated 
spot areas also follow a log-normal distribution. This finding was 
confirmed by later studies (e.g. Baumann & Solanki 2005; Hath 


away & Choudhary 2008] |Kiess et al.||2014| >. Moreover, 


Bau- 


mann & Solanki 1 2005|> have shown that the log-normal distri¬ 


bution also nicely fits the observed distributions for the instan¬ 
taneous area and for the maximum area of sunspots. However, 
according to Jiang e t al.| ( j20TI| , the area distribution for groups 
smaller than 60 MSH (millionth of the solar hemisphere) is bet¬ 
ter described by a power law. 

During its life, a given sunspot group grows until it reaches 
a maximum area and then decays. The growth ( V F) and decay 
(T) rates, i.e. the time derivative of the group’s area during each 
of these phases, are found to be dependent on the group areas, 
the activity cycle, the phase of the cycle, and the latitude (e.g. 
Moreno-Insertis & Vazquez| 1988| Howard 1992| Petrovay & 
van Driel-Gesztelyi|1997||Hathaway & Choudhaiy 2008) Javara- 


iah|2012 i. Small groups grow faster than they decay, while the 

growth rates of large groups are smaller than their decay rates 
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(Howard] | 1992| . |Hathaway & Choudhary | ( j2008] > found a linear 
relation between the decay rates and the group area, but the ero¬ 


sion model proposed earlier by Petrovay & van Driel-Gesztelyi 
(1997) indicates a non-linear relation of the type T oc A u ;l , where 
A is the group area at a given time. More recently, Javaraiah 


(2012)1 suggested that the relation between the decay rates and 


areas could be better described by a power law of the form 
T = exp(yi )A n . The constant j 2 was found to vary from ~ 0.45 
to ~ 0.70, when considering individual cycles and different 
phases of each cycle (with A expressed in MSH). When assum¬ 
ing the sunspot data from 1874-2011, Javaraiah (2012 1 found 
that yi ~ 0.26 and yi ~ 0.613. On the other hand, the same 
study did not produce conclusive results regarding the relation 
between the group’s growth rate and its area. 

The areas of sunspots and sunspot groups are also related 
to their lifetimes, which can range from hours to months, de¬ 
pending on their size. This dependency is described by the 
Gnevyshev-Waldmeier (GW) rule (Gnevyshev 1938 1 Waldmeier 
1955)1, according to which, 


A-m - D GW T . 


( 2 ) 


1938) Filing & Schwentek 1992[|Sabarinath & Anilkumar|2008f 

Volobuev||2009[ Du |2011) i. In particular, |Hathaway, Wilson, &] 

Reichmann| ( 1994| ifound that the observed number of sunspots 
is nicely fitted by the following function of time, 


Here, A m is the sunspot or sunspot group maximum area (in 
MSH), T is the corresponding lifetime, and is a constant 
of proportionality (around 10 MSH day '). The determination of 
a precise value for Dow is hampered by the difficulty in mea¬ 
suring the spots’ lifetimes due to the nightfall, the solar rotation 
(lack of observations of the invisible side of the Sun), and limb 
darkening (e.g. |Henwood et ak 2009[|Blanter et al.|2006)|Solanki| 
2003) 1. In spite of these difficulties, some studies have been car¬ 
ried out, indicating that Dqw might be larger than first esti¬ 
mated. Petrovay & van D r i c I - G e s z t e I y i j( l 997) found that Dow is 
10.89 + 0.18 MSH day -1 for individual sunspots. More recently, 
Henwood et al. ( 2009) studied long-lived sunspot groups and es¬ 
timated that Dgw = 11.73 + 0.26 MSH day" 1 . 

As was first reported by Carrington (18631, the sunspot for¬ 
mation latitude also varies periodically with time. At the begin¬ 
ning of a new cycle the first spots appear at latitudes of about 
±40°. The succeeding spots form at progressively lower lati¬ 
tudes, being rarely observed within ±5°. At the solar minimum 
the last spots of the cycle emerge at low latitudes, while spots of 
the new cycle start to form at high latitudes. This behaviour is 
known as the Sporer law and it may be seen in the butterfly di¬ 
agram (or Maunder diagram; |Maunder| 1904j ). Despite the pres¬ 
ence of short plateaus at intermediate latitudes (~ 10° - around 
maximum; |Chang||2012 1 , |Hathaway] ( 12011) 1 found that the drift 
of the sunspot zones follows an exponential function, where the 
average latitude, L$, is given by 


L s (t) = Lo exp 




(3) 


where Lo is the mean latitude at the time to, and t is expressed 
in years. By considering the intermediate phases of the cycle 
where there is no overlap between consecutive cycles, Jiang et al. 
( |2011] > verified that the evolution of the average latitude can also 
be described by a second-order polynomial. 

Chang ( |2012) > found that the spatial distribution of the 
sunspot groups at each time t is bimodal and that it can be 
described by a double Gaussian; instead, Ivanov et ak] ( |2011) > 
showed that one single Gaussian describes the data reasonably 
well. Moreover, the width of the sunspot formation zone, cr L , 
also varies over the solar cycle (Gleissberg 1958). According to 
Miletskii & Ivanov] ( |2009| ), |Ivanov et al.| ( |20lT| ] and |Ivanov &] 
Miletskii (201 1)>, this width is a linear function of the activity 
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level. However, according to Jiang et al. ( |201 1] >, it can be de¬ 
scribed in relation to the average sunspot group latitude by a 
second-order polynomial of the form 


ox 

Ls 


+ C, t 



(4) 


where a (r , b ir , and c rr are the coefficients of the polynomial, t mln 
corresponds to the minimum, and P c is the period of the cycle. 

We note that sunspots are depressed (Wilson depression), 
thus they move according to the subphotospheric layers, i.e. 


slightly faster than the solar surface (e.g. Zappala & Zuccarello 

1991 

Zuccarello 1993( 

Abuzeid & Marik 1997 ; Schou et al. 

1998 

Kitchatinov 2011 

). Since the spots’ depths decrease as 


they evolve, younger spots also move faster than the older ones. 

Finally, we note that the properties described above vary 
from cycle to cycle and some are found to be correlated (e.g. 
Solanki et al.|2002 |2008 1 [Hathaway |20 1 0| Jiang et al.|201 1 1. 


3. Empirical solar cycle model 

The primary goal of this work is to produce a tool capable of 
reproducing an activity cycle that retains the main observational 
properties of the solar cycle. To that end, we develop an empir¬ 
ical model to generate sunspot groups as a function of time and 
gradually adapt the model assumptions and inputs until our goal 
is reached. To decide whether or not a given assumption/input is 
a better representation of the observational data than the previous 
ones we use the Kolmogorov-Smirnov test, which is described in 
section l4~2l 


3. 1. Number of sunspot groups 


In our model, each sunspot group is generated independently 
from the others. This could not be assumed if we were consider¬ 
ing individual spots, since the formation of spots within the same 
group is not independent. 

At each time step (fixed on one day to be comparable to the 
daily records of the sunspot data), N groups are formed. The 
number of generated groups is randomly determined using a 
Poisson distribution with a mean value N m that depends on time. 
In the current version of the model, N m is taken to be one sixth 
of /V s , where the function of time A' s is derived from a fit of equa¬ 
tion o to the number of observed sunspot groups for solar cycle 
23 (Fig .0. With this choice for ;V rn , we find that the function N s 
derived from the fit to the synthetic data is always in reasonable 
agreement with that derived from the solar data. 

Other ways to determine N s from the solar data were ex¬ 
plored. The functions used by Du ( 201 1[ > and Sabarinath & 
Anilkumar (2008) (the latter takes into account the double peak 
feature of the solar cycle) led to results that are comparable to 
those obtained with equation <0- 


3.2. Sunspot formation zone 

In our model, the latitude of each spot, L, is determined ran¬ 
domly through a single-Gaussian distribution, with a mean lati¬ 
tude Ls and dispersion <r L , both dependent on time. For Ls, we 
use Hathaway’s exponential function given by equation (|3]i fit¬ 
ted to the northern hemisphere solar data for cycle 23 (Fig. [2}, 
with to fixed at the value found in section |3.1| We opted to use 
solar data from a single hemisphere because in our model we do 
not account for the north-south asymmetry. For cr L , we assume 
a second-order polynomial (Fig. |3j, as suggested by Jiang et al. 



Fig. 1. Variation of the number of sunspot groups over the solar cycle 
23 (black). The red smooth line corresponds to the fit obtained using 
the function defined by Hathaway et al.| ( p994} >. The observational data 
were extracted from NOAA. 


( |2011| >. Other functions were considered to describe the latitu¬ 
dinal distribution of the sunspot groups (both for Ls and cr L ). 
However, those led to a lower level of agreement between the 
synthetic and the observational data. 


~o 

zs 


a 



1998 2000 2002 2004 2006 

Time (years} 


Fig. 2. Comparison between the latitudinal distribution of sunspots and 
the exponential fit of |Hathaway| l |2011] ). The observational data was ex¬ 
tracted from NOAA. 



1998 2000 2002 0004 2006 
Time (years) 

Fig. 3. Ratio between the standard deviation of the latitudinal distri¬ 
bution and the corresponding mean latitude (observational data from 
NOAA). The red solid line corresponds to the best fit. 
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3.3. Sunspot group areas and lifetimes 

Having defined the procedure to generate the number of sunspot 
groups and their position in latitude at each time step, we need 
to associate a maximum area ( A m ) to each of them. Although 


the mean group area varies over the solar cycle (e.g. Jiang et al. 
|2011] >, we fix the area distribution. In accordance with the dis¬ 
cussion in section [2] in our model the sunspot group maximum 
areas are drawn from a log-normal distribution whose param¬ 
eters are obtained from a log-normal fit to the sunspot groups 
observations for cycle 23, considering each group only once and 



0 100 200 300 400 500 600 

Group area (MSH) 


Fig. 4. Accumulated (blue) and maximum (green) area distributions 
for the solar data from NOAA. The black lines are the respective log¬ 
normal fits. 


With the sunspot maximum area in hand, we can in principle 
determine the group’s lifetime through the GW rule. However, 
we have found that when the GW rule is taken for all ranges of 
areas, the accumulated area distribution obtained with the em¬ 
pirical tool is not in agreement with the observed distribution: 
the number of small groups in the synthetic distribution is lower 
than that found in the Sun and the peaks of the two distributions 
do not coincide. The upper panel of Fig. [5] illustrates this dis¬ 
agreement. 

The GW rule is based on the observation of long-lived spot 
groups (e.g. groups that live longer than ~ 20 days; |Henwood| 
|et al.||2009] l. Hence, it is possible that this rule may not be ad¬ 
equate for the smallest groups. In fact, by tracking the small 
groups in the solar data for cycle 23 and comparing the time they 
remain visible with the lifetime predicted by the GW rule, it is 
possible to verify that the GW rule generally underestimates the 
lifetimes of the small groups. Moreover, there is a significant dis¬ 
persion around the area-lifetime relation, which for the smallest 
groups is strongly asymmetric (since lifetimes cannot be neg¬ 
ative). With this in mind, we have checked whether increasing 
the lifetimes for the small groups would improve the agreement 
between the observed and synthetic accumulated area distribu¬ 
tions and found that substituting the GW linear relation by an 
exponential relation at the lower areas end, our model produces 
an accumulated area distribution that is in better agreement with 
the observed one (lower panel of Fig. [5]>. The modification of the 
GW rule for groups with areas smaller than 85 MSH used in our 
model is 

T = 5 exp(6.2591 x l(T 3 A ffl ). (5) 

Figure[6]shows the current area-lifetime relation used in our em¬ 
pirical model (red line). 

1 The bin for the smallest groups is incomplete, as the sunspot records 

are limited to areas larger than 10 MSH. 
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Fig. 5. Accumulated area distribution for synthetic sunspot groups when 
assuming the GW rule for all areas (upper panel) and when using the 
modified area-lifetime relation (lower panel). The black, solid line is 
the log-normal fit to the observational data.F] 



Fig. 6. The red line shows the area-lifetime relation assumed in the cur¬ 
rent model while the blue dashed line shows the lifetime predicted for 
small groups by the GW rule. 


Taking the group’s maximum area and lifetime, we deter¬ 
mine the group’s area at a posterior time (while A > 0) by ap¬ 
plying a decay rate of the form T = exp (y\ )A yi . We then assume 
that the period during which the group grows corresponds to the 
difference between its lifetime and the decay time. In the absence 
of a well-established relation between the group’s growth rate 
and its area, we opted for a power law to describe that relation, 
i.e. T = exp (if/] )A' //2 . We started by considering the values found 
by Javaraiah ( 2012} ) for j\ and y 2 , making y\ = d\ = 0.26 and 
y 2 = if /2 = 0.613. These parameters were progressively changed 
until a reasonable agreement between the observations and the 
synthetic data was reached. The current version of the model 
considers T 1 * = exp (0.17) A 0 46 and T = exp (0.17) A 0 47 , which is 
consistent with the fact that large groups have higher decay rates 
than growth rates. For small groups the growth and decay rates 
given by the expressions above are essentially the same. How- 
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ever, since the growth time is taken as the difference between the 
lifetime and the decay time, in practice the smallest groups do 
not show a growth phase in the daily records, which might be 
interpreted as a fast growth, where they reach the maximum area 
in a time shorter than the interval between consecutive records 


(one day). We also tested the linear relation found by Hathaway 
|& Choudhary| ( |2008| , but we found that decay rates are too high 
when compared with the lifetime of the group. Assuming the 
linear relation for the growth and decay rates, the time interval 
from the first appearance (with A ~ 0) to the last appearance 
(with A ~ 0) is much shorter than the lifetime from the modified 
GW rule. 


3.4. Sunspot visibility 

To reproduce the daily sunspot records of the Sun, we take into 
account the solar rotation and the fact that spots are observed 
only when they are on the Sun’s visible side. The group’s lon¬ 
gitude is determined randomly from a uniform distribution be¬ 
tween 0 and 2n. If this quantity is smaller than n, we consider 
that the sunspot group is on the visible side, otherwise we con¬ 
sider that it cannot be observed. Taking into account that the ro¬ 
tation velocity of the subphotospheric layers is not very different 
from that of the solar surface, we assume the following param¬ 
eterization of the groups’ rotation velocity, to, as a function of 
their latitude, L ( |Snodgrass|1983||Snodgrass & Ulrich|1990| > 

to(L) = 14.71 -2.33 sin 2 L- 1.78 sin 4 L. (6) 

Spot groups that emerge on the visible side of the Sun can move 
towards the invisible side and then eventually become visible 
again depending on their lifetimes and on the solar angular 
velocity at the latitude they emerge. Moreover, groups that 
emerge on the invisible side of the Sun can become visible. Both 
these facts are taken into account in our model and only groups 
with an area greater than ~ 10 MSH are considered visible (in 
analogy to the sunspot data). 

A schematic summary of our empirical model is shown in 
Appendix [A| 


4. Results 

4.1. Synthetic data 

The synthetic data produced with our empirical model provide 
information about each generated group. In analogy to what is 
done in the NOAA databases, the code yields the sunspot group 
records, which include date, latitude, group area, lifetime, posi¬ 
tion in longitude, rotation rate of the solar surface at the group’s 
latitude, and an identification number. These records are then 
used to compare our results with observed data. 

A comparison between the number of sunspot groups ob¬ 
served over cycle 23 (in black) and those obtained in one real¬ 
ization of our model (in red) is shown in the left panel of Fig. [7] 
The similarity in the shape and spread of the two curves is quite 
evident. The same similarity is found when comparing the total 
group areas, i.e. the total area covered by sunspot groups in each 
day, (Fig. □ right panel) and the group latitudes (Fig. [8]) for the 
real and synthetic data. 

4.2. Comparison test 

In order to test and improve our model we have quantified how 
closely related the observed and synthetic data sets are by apply- 



-60 


1998 2000 2002 2001 2006 

Time (years) 

Fig. 8. Same as in Fig. [ 7 ] but for the real and synthetic latitudinal distri¬ 
bution of sunspot groups. 


ing the Kolmogorov-Smirnov test ( Kolmogorov| 1933 Smirnov 
1939) . This test compares the cumulative distribution functions 
of two samples, using the maximum deviation between them. 


O n| .n 2 = max \F\ Al {x) - F x n 2 (*)|, 


(7) 


where n \ and no are the number of elements of sample 1 and 
sample 2, respectively, and F t and F 2 are the corresponding cu¬ 
mulative distribution functions. 

The null hypothesis - i.e. that both samples result from the 
same distribution - is rejected at significance level a if 


D„ 


> c(a ) 


n 1 + n 2 
n 1 n 2 


( 8 ) 


where c(a ) is a constant that depends of the significance level to 
be considered. 

More than to reject the null hypothesis at a given level a, in 
our case this test was used to identify the aspects of the model 
which needed to be improved. The synthetic sunspot cycles ob¬ 
tained from early versions of the model resulted in large values 
of Z) ni indicating that they did not provide a good descrip¬ 
tion of the observed properties of the sunspot cycle. By compar¬ 
ing the observed and synthetic cumulative distribution functions 
for the total area covered by sunspot groups and for the groups’ 
latitudinal distribution, we could decide where and how to im¬ 
prove our model. An example of this is provided in Fig. [9] where 
we compare the cumulative distribution functions for the total 
group areas of the synthetic and real data. In this case, F cor¬ 
responds to the fraction of days with total area below a given 
value, and n\ = n 2 = 4017 is the total number of days consid¬ 
ered in the real and synthetic data. The left panel corresponds to 
a sunspot cycle reconstruction obtained with the current version 
of the model that considers an exponential function (equation 
(|5]i) to correct the GW rule for groups smaller than 85 MSH, 
while the right panel shows the results for a reconstruction ob¬ 
tained when adopting the GW rule for all group areas. Small 
groups from the former reconstruction live longer than groups 
with similar areas from the latter. This leads to an increase in the 
daily number of sunspot groups and, consequently, to a larger 
total area covered by the groups in the reconstruction obtained 
with the model that incorporates the corrected area-lifetime re¬ 
lation than with the other. The result is a shift of the cumulative 
distribution function towards larger areas and a better agreement 
with the observations. 
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Fig. 7. Comparison between real solar data for cycle 23 (black) and synthetic data (red). Shown are the number of sunspot groups (left panel) and 
the total group area (right panel). 




Fig. 9. Cumulative distribution functions of the synthetic (red) and real data (black) for the total group areas of a reconstruction obtained from the 
current version of the model (left panel; £) ni , n , ~ 0.0141) and a reconstruction obtained when assuming the GW rule for the area-lifetime relation 
(right panel; D ll|n2 ~ 0.0798). 


Figure[K)]compares the cumulative functions for the real and 
synthetic latitudes resulting from the same reconstruction, ob¬ 
tained with the current version of the empirical model. Here, 
the cumulative distribution functions, F, represent the fraction 
of groups that become visible (first appearance) at a latitude 
lower than a given value, and n t = 2801 and n 2 = 2919 cor¬ 
respond to the total number of different observed groups in the 
real and synthetic data, respectively. Despite the small value of 
the statistics Z) ni n ,, the difference between the two distributions 
is significant and its interpretation is relatively straightforward: 
the cumulative distribution function for the solar data indicates 
that the southern hemisphere retains almost 55% of the visi¬ 
ble sunspot groups, while for the synthetic data the groups are 
more evenly distributed by the two hemispheres. This discrep¬ 
ancy results from the hemispheric asymmetry that is known to 
be present in the data, but that is not accounted for in our model. 
To verify this, we compared the cumulative distribution func¬ 
tions for the absolute values of the latitude (Fig. m finding a 
lower value of the statistics D„ u „ 2 . 

The results from the KS-test for the synthetic sunspot cycle 
discussed above are summarized in Table |T] At the significance 
level a = 0.1 the null hypothesis for the total group area and ab¬ 
solute group latitudes is not rejected. Although the non-rejection 
of the null hypothesis does not allow us to conclude about its 
veracity, it certainly reinforces the expectation born from the 
direct inspection of Figs. [7] and [8] that the synthetic cycles ob- 
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Fig. 10. Cumulative distribution functions of the synthetic (from the 
current model; red) and real data (black) for the spot latitudes. The max¬ 
imum difference is D ni ,n 2 ~ 0.0377 (blue). 

tained from our model retain the main observed properties of the 
sunspot cycle. 

As the results from our empirical model are stochastic, we 
can perform Monte Carlo simulations to obtain the distribu¬ 
tions for the statistics While n 2 is constant for the to¬ 

tal area covered by sunspot groups, for the latitudinal distribu¬ 
tion ii 2 varies from reconstruction to reconstruction and accord- 
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n\ 

«2 

^n\,ri2 

D(a = 0.1) 

Total group area: Fig. 

9 left 

4017 

4017 

0.0141 

0.0272 

Latitudes - FTg. 
Absolute Latitudes - Fig. 

10 

11 


2801 

2919 

0.0377 

0.0283 

0.0323 


Table 1. Results from the KS-test for the synthetic sunspot cycle discussed in the text. The ri\ and n 2 are the sample sizes for the real and synthetic 
data, respectively. D = c(0.1) V(”i + 'ril/Oqui) is the right-hand side of equation |8]l for a significance a = 0.1. 



Fig. 11. Cumulative distribution functions of the synthetic (from the 
current model; red) and real data (black) for the absolute values of the 
group latitudes. The maximum deviation is -D ni ,n 2 ~ 0.0283. 



1 2 3 '1 5 

C „,, n2 — Total group- area 


Fig. 12. Distribution of C niA2 = D niiIl2 / VOq + n 2 )l(nin 2 ) for the total 
area covered by sunspot groups, when assuming the GW rule for all 
areas (blue histogram) and when using the modified area-lifetime rela¬ 
tion (red histogram). The vertical lines indicate the levels of significance 
a = 0.1 (dotted), a = 0.05 (dashed) and a = 0.001 (dash-dotted). 


ing to the case considered (both, southern or northern hemi¬ 
spheres). Thus, rather than considering the distributions for 
D nu n 2 , we consider those for the « 2 -independent quantity C ni , n2 = 


12 


and 


13 


summarize the re- 


Ai,,n 2 / V(«t +«2)/(»i«2)- Figures _ __ 

suits obtained from 5000 cycle reconstructions. The distribution 
of C nii n 2 obtained from the analysis of the total group area when 
the GW rule is assumed for all ranges of area (Fig. [12] blue 
histogram) is shifted towards larger values of C ni ,n 2 than that 
obtained from reconstructions that apply the correction to the 
area-lifetime relation for small groups (Fig. [T2] red histogram). 
This confirms that the modified area-lifetime relation used in 
our model produces results that are statistically in better agree¬ 
ment with the solar data. Concerning the sunspot groups’ lat¬ 


itudes (Fig. 13 i, the comparison of the distributions for C n 


clearly confirms our earlier findings. The consequence of the 
non-inclusion of the hemisphere asymmetry in our model is that 
our synthetic cycles compare significantly better when consider¬ 
ing the absolute values of the latitude. 


The properties of the sunspot emergence are cycle- 
dependent, stronger cycles usually having higher latitudes and 
wider sunspot formation zones than weaker cycles (e.g. Solanki 
et al. 2008| Jiang et ak]|2011 1 . With this in mind we have ap¬ 
plied our model to a second cycle, namely cycle 22, which is 
stronger and more asymmetric than cycle 23. In agreement with 
the works mentioned above, we found that the average latitude 
and the width of the sunspot formation zone obtained in our re¬ 
constructions for cycle 22 are larger than those obtained for cy¬ 
cle 23. Although the results from the KS-test for the latitudes 
were found to be worse for cycle 22 owing to the hemispheric 
asymmetry, we found that with regard to the total area covered 
by spots and to the absolute latitudes the synthetic data for this 
cycle is also in qualitative agreement with the real sunspot data. 



1.0 1.5 2.0 2.5 3.0 

C„ 1-n , - Latitude 


Fig. 13. Distribution of C ni ,„ 2 = D aiA2 / VOri + >h)/( n i n 2 ) f° r the lat¬ 
itudinal distribution: from the group latitudes (blue) and absolute lat¬ 
itudes (red). The vertical lines indicate the levels of significance a = 
O.l(dotted), a = 0.05 (dashed) and a = 0.001 (dash-dotted). 


5. Conclusions 


In this work, we have presented an empirical tool for the stochas¬ 
tic reconstruction of sunspot cycles. With the parameters adopted 
in the version presented here, our tool produces synthetic daily 
sunspot records which retain the main properties of the real solar 
data. 

A crucial assumption of our model is that different sunspot 
groups are generated independently. Despite evidence that 
sunspots tend to form within active longitudes (e.g. |Bumba &| 
Howard|1965[ Bogart|1982| |Jiang et al.|201 l| i, pointing to pos¬ 
sible correlations between their emergence, we found that the 
properties of the sunspot cycle are reasonably well reproduced 
under this model assumption. If a significant dependence be¬ 
tween the generation existed we would expect that to have been 
reflected, for example, in a lack of agreement between the dis¬ 
persion of the number of synthetic and observed sunspot groups 
as a function of time. In fact, we have checked that this problem 
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would occur if individual spots were (incorrectly) considered as 
independent events. 

An important by-product of our work was the verification 
that the GW rule is only appropriate for the largest sunspot 
groups. For the smallest sunspot groups, this rule seems to un¬ 
derestimate the groups’ lifetime. We have proposed a modified 
area-lifetime relation for small groups which leads to a closer 
agreement of the synthetic sunspot cycle with observations. 

The quantity, quality, and diversity of solar data, can only 
be adequately reproduced by a relatively complex model that in¬ 
cludes a number of empirical parameters and functions describ¬ 
ing the average observed properties of the sunspot number, area, 
latitude and rotation. In contrast, for other stars the observational 
constraints are much more limited. Thus, the application of our 
tool to the study of activity cycles on other stars will require 
the identification of the model parameters that have a significant 
impact on the activity-related stellar observables, including the 
frequency shifts (work in progress). Another potential applica¬ 
tion of this kind of tool is related to the search for exoplanets, 
where new strategies for reducing activity signatures in the ra¬ 
dial velocity and transit observations can be designed. 
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Appendix A: Schematic overview of the empirical model for the solar cycle 

Figure |AJ~| provides a schematic summary of the procedure underlying our empirical tool. The green box illustrates the detailed 
treatment of the evolution of a sunspot group. All the parameters included in the model have been found to be important in order to 
obtain results that are in reasonable agreement with the observations. In the following scheme, U, Pois, AT, and In AT are, respectively, 
the uniform, Poisson, Gaussian, and log-normal distributions. The location in longitude of each sunspot group is represented by O. 


/ ' 

P c = cycle period 

dh = step in time 
tj = initial time 
t m in = I s ' mininum 


Temporal evolution (t-years) 




1 .N x = 


Qi(t-tp) 

(f— t 0 ) 2 lb 2 , 

e c ± 

-(t—f 0 )/7.5 


, t-t min \ t~t min 
3. a,=a,.+b,. - +c_ - 



Fig. A.l. Summary of the empirical cycle model. The green box shows how the evolution of a spot group is determined in our model. 
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